{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bb07c35f",
   "metadata": {},
   "source": [
    "## 周期性成分的显着性检验"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1c18518",
   "metadata": {},
   "source": [
    "   Fisher 的g统计量评估白噪声中正弦分量的重要性。Fisher 的g统计量是最大周期图值与 1/2 频率间隔 (0, Fs/2) 上所有周期图值之和的比率。g统计量和精确分布的详细描述可以在参考文献中找到。\n",
    "   在高斯白噪声中创建一个由 100 Hz 正弦波组成的信号，均值为零，方差为 1。正弦波的幅度为 0.25。采样率为 1 kHz。将随机数生成器设置为可重现结果的默认设置。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4117f7a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fs = 1000\n",
    "t = np.arange(0,1,1/fs)\n",
    "rng = np.random.default_rng()\n",
    "x = 0.25*np.cos(2*np.pi*100*t) + (2*rng.random(size=t.shape)-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e80298eb",
   "metadata": {},
   "source": [
    "使用获得信号的周期图periodogram。排除 0 和奈奎斯特频率 ( Fs/2)。绘制周期图。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ecce8b5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "F, pxx = signal.periodogram(x, fs,window='boxcar',nfft=len(x))\n",
    "# python boxcar等于matlab中retwin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "31d4705e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB2cklEQVR4nO2deZwcRfn/P0/P7JE9cockJCEJEEg4wpUEkGuR+xLx9isKoiCKfvXngSDqFxUURQFPEC8uFVARUe5ANiFCyEkSEhLIRe472WTvnZn6/VFV3dXV1T09s3Nsduv9eu1rp+/q7up66jnqKWKMwWKxWCwWFafcBbBYLBZLz8MKB4vFYrEEsMLBYrFYLAGscLBYLBZLACscLBaLxRLACgeLxWKxBLDCoQdARM8S0VV5HruOiM4tdJksByZENI6IGBEly10WACCiq4lotrLcTESHlrNMOkR0BhGtLHc5ehpWOOSJaJTbRGXfRkQPEFFdPudijF3EGHuw0GXsyRDRrUTUJZ7fXiJ6lYhOLXe5JERUSUQ/I6KNoozriOieIl/T15AW6Rqy3u5Xnvv1RFSStoAxVscYWyPK8gAR3ZbvubRvUP79KsZxjIgOV8r0CmPsyHzLkeVa3brHcmKFQ/e4jDFWB+BEAFMAfDuXg4nTI98BESVKcJnHxPMbBmA2gCeIiEpwXR8hveybwd/pNAD1ABoALCxhsYwU6L1cxhirBzAWwB0AvgngDwU4bzm4TAgc+ffFcheo18AYs395/AFYB+BcZflOAP8Rv08B8CqAvQAWA2hQ9msEcDuA/wJoA3C4WPdZsd0BFzLvAtgO4CEAA5TjPym27QJwi1oOAFUA7gGwWfzdA6BKOfZGAFvEts8CYAAOF9seAHAvgGcAtAA4F8AlABYB2AdgA4BblXONE8d/WmzbA+B6AFMBLBH3/quI53crgEeU5aPF+YYCOBjAUwB2A1gF4FqxT7V4ZkPF8i0AUgD6i+UfALhHeRY/BbAewDYA9wHoJ7Y1ANgI3ihuBfCwoXz/AfCVLO//ZgDLxb3/CUC1sv1SAG+I5/AqgMnKtjEAngCwQ7zHXwGYBKAdQBpAM4C93XwvyTj1VqybBiAD4Jgcnt3XwOvnFgCfVs41RLy7fQDmincyW9nOwOv8dQC6AHSK+/03gG8A+IdWtl8A+Hnce1G2HQ5gJoAmADvBOyIAMEuUoUVc96PynrTzfgO8HreAC87hAJ4FsB/AdACDlP3/JupRkzj/0WJ94B7F+oMB/EO8/7UA/rfc7ZnxGZa7AAfqH/yN8hgAy8SHMEp88BeDN/TnieVhYt9G8dEdDSAJoAJ+4XANeIN4KIA68EbkYbHtKFHJzhQf8F3gjaMsx/cBzAFwEHhv/FUAPxDbLhQV+GgANQAeQVA4NAE4TZS7Wnw0x4rlyeANxfvF/uPE8feJfc8Hb9yeFNcfBd54nBXy/G6FEA7iXu4EsF4szwLwG3He48VH9F5l2wfF7xcArAZwkbLtCvH7bvBGajB4z//fAH4ktjWI5/Zjce1+hvJ9W7ynL4hnQIb3/6Z494PBhf1tYtsJ4t5PBpAAcJXYv0osLxblqxX3eLo47mooDWk330ts4SDWrwfw+Rye3ffB6+7FAFohGksAjwJ4XNzbMQA2wSAclHu7Tdk2ErwxHiiWk+I5npTLvYhtfwXvPMhndrqpDMo96cJhDrhAkPV4oXiv1QBeBvB/yv7XiOckO2dvaO9PvUcHwAIA3wVQCf6drwFwQbnbtMAzLHcBDtQ/UYGawXuG74I3Zv3Ae6MPa/s+D+Aq8bsRwPe17Y3whMNLAL6gbDsSvPeRFBXqUWVbLXivRAqH1QAuVrZfAGCd+P1H+YGL5cMNH+pDWe75HgB3i9/jxPGjlO27AHxUWf4HQnrf4MKhUzy/7eKDOwm8sU0DqFf2/RGAB8TvH4D3JpPgwu7L4KYRqVUMAUDgjcxhyjlOBbBW/G4Q166OuNcEgBvAG/0OcG3rKu39X68sXwxgtfh9L4RQVravBHCWKMcOGBpvhAuHfN5LrsJhDnhjGufZtannF+/vFPHMugBMVLb9EDGFg1j3LDxN8VIAy2N+g/JPHvsQgPsBjDYcF0c4fEKrx/cqy18C8GRImQaK8w8w3SN4h2G9dszNAP4U9Y7L8dcj7d0HEO9njA1kjI1ljH2BMdYGbsf9sHD27SWivQBOB+8VSTZEnPNgcGEjeRe8IRwutrnHMsZawBvkqGMPVrap1zWVwbeOiE4mohlEtIOImsDNRkO1Y7Ypv9sMy1FO+sfF8zuIMfZextgCUc7djLH92n2MEr9ngn/MJwJYCuBF8Eb3FACrGGO7wLWmGgALlHfwnFgv2cEYaw8rGGMszRj7NWPsNPAP/nYAfySiScpu6vNSn/VYAF/T6sAYsX0MgHcZY6mI56KTz3vJlVHgZrw4z26XVv5W8Pc8DLyu6s8lFx4EcKX4fSWAh7PsL79B+fc7sf5GcEE3l4iWEdE1OZYjVr0mogQR3UFEq4loH7hgAcLfx1gAB2t141vg33ePwgqHwrMBXHNQK2wtY+wOZR8Wcfxm8AokOQRcjd8Gbt8dIzcQUQ14Tznq2M3i9xYAo5VtYxBEL9dfwM0LYxhjA8BNSMV2GG8GMJiI6pV1h4CbJwBuKjsSwBUAZjLGlovtF4MLDoDbmNvAbb/yHQxg3PktiXoHPhhjbYyxX4P7Fo5SNqnPUH3WGwDcrtWBGsbYX8W2Q0Kc4GFlKup7IaKp4MJhNuI9uzB2gNdV/bmEYbrfJwFMJqJjwDWHP8e4bvDEjG1ljF3LGDsYwOcA/EaNUCog/wPgcnBf0ABwzQ3w3od+jxvAtTC1btQzxi4uQtm6hRUOhecRAJcR0QWiV1FNRA1ENDrrkZy/Avh/RDRehMb+ENyZlgLwdwCXEtHpRFQJbvd1tGO/TUTDiGgouBnqEbHtcQCfJqJJQqh8J0ZZ6sF78e1ENA38QygqjLEN4ALgR+LZTQbwGYj7YIy1gttsb4AnDF4F7z3PFPtkAPwOwN1EdBAAENEoIrogbjmI6CvivfUjoqQYh1IP7giW3EBEo4loMLhJ5jGx/ncArhc9fCKiWiK6RAi8ueCC+g6xvpqIThPHbQMwWrzbKAryXoioPxFdCu4neIQxtrQ7z44xlgb3kd1KRDVEdBS4vyWMbeA2d/Uc7eD1/C8A5jLG1ud5bx9Wvrk94I10Juy63aAe3Oy4C1zj+qG2Xb/WXAD7ieibom4liOgYIaB7FFY4FBjRuF0OriruAO8pfAPxn/UfwVXpWeCRDO3gNk4wxpaBN4p/AW9g9oBHjkhuAzAfPMpiKbgT7TZx7LPgtvoZ4A7vOeKYjoiyfAHA94loP7igeTzmPXSXj4P3wDYD+Ce482+6sn0muDN0rrJcD/7MJN+EuE+h7k8H1zji0grgZ+B+jZ3gz/2DTMToC/4C7hRfA+7vkc96PoBrwaOQ9ohyXC22pQFcBu7zWQ/+/j4qzvcyeGDDViLaGVG27r6Xf4tjN4ALtbvAo84k3Xl2XwQ3uWwFt7f/KWLfPwA4SphXnlTWPwjucM9mUgL4vajjHP4p1k8F8DoRNYNrWV9W3t2tAB4U1/1IzPsK4yFw09km8Mi1Odp23z2K938peKDFWvC69XtwraNHQcIhYuljCNv5m+ChrrnYvy3gA7DAgwimZ9vXkhtEdAiAFQBGMMb2lbs8fRWrOfQhiOgKIqoiokHgYZz/toLB0pMQg0K/Ch6VZwVDGbHCoW/xOfCww9Xg4aKfL29xLBYPIqoFHzx3HoD/K3Nx+jzWrGSxWCyWAFZzsFgsFkuAHpHWt7sMHTqUjRs3Lq9jW1paUFtbW9gC9XDsPfcN7D33DbpzzwsWLNjJGBtm2tYrhMO4ceMwf/78vI5tbGxEQ0NDYQvUw7H33Dew99w36M49E1HoCHZrVrJYLBZLACscLBaLxRLACgeLxWKxBLDCwWKxWCwBrHCwWCwWSwArHCwWi8USwAoHi8VisQSwwqGP0bhyO3a2ZbLvaLFY+jS9YhCcJT5X/2keqhLAhy4qd0ksFktPxmoOfZCOdLlLYLFYejpWOFgsFoslgBUOFovFYglghYPFYrFYAljhYLFYLJYAVjj0IeysfxaLJS5WOPQhMlY2WCyWmFjh0IfIWM3BYrHExAqHPkTaqg4WiyUmVjj0IaziYLFY4mKFQx/CmpUsFktcrHDoQ1jhYLFY4mKFQx8iY5OxWiyWmFjh0IewmoPFYomLFQ59CCscLBZLXKxw6EPYSFaLxRKXsggHIrqTiFYQ0RIi+icRDVS23UxEq4hoJRFdUI7y9Vas5mCxWOJSLs3hRQDHMMYmA3gbwM0AQERHAfgYgKMBXAjgN0SUKFMZex1WOFgslriURTgwxl5gjKXE4hwAo8XvywE8yhjrYIytBbAKwLRylLE3Ys1KFoslLj1hDulrADwmfo8CFxaSjWJdACK6DsB1ADB8+HA0NjbmdfHm5ua8jz3Q2NHqxbL2lXuW9KX3LLH33Dco1j0XTTgQ0XQAIwybbmGM/UvscwuAFIA/53p+xtj9AO4HgClTprCGhoa8ytnY2Ih8jz3QeHdXCzCrEQD6zD1L+tJ7lth77hsU656LJhwYY+dGbSeiqwFcCuAc5k00sAnAGGW30WKdpQBYs5LFYolLuaKVLgRwI4D3McZalU1PAfgYEVUR0XgAEwDMLUcZeyPWIW2xWOJSLp/DrwBUAXiRiABgDmPsesbYMiJ6HMBycHPTDYyxdJnK2OuwM8FZLJa4lEU4MMYOj9h2O4DbS1icPkPa5layWCwxsSOk+xDWrGSxWOJihUMfwgoHi8USFysc+hBWNlgslrhY4dCHsHNIWyyWuFjh0IewZiWLxRIXKxwKRDrDsGN/R7mLEYlVHCwWS1yscCgQP35uBabePh27WzrLXZRQrOZgsVjiYoVDgXjuza0AgH1tXWUuSTgZqzpYLJaYWOFQIGSvPOFQmUsSjpUNFoslLlY4FAjZK6eeKxts+owikckwdKbs8HNL78IKhwKREsIh6fTcR5q2wqEofOnRRTji28+WuxgWS0HpuS3ZAYY0K/VkzcGalYrD00u2lLsIFkvBscKhQMiGtyd3zm20ksViiYsVDgVCjj7uyQ2w9TkUFxsNZulNWOFQIKRw6MnNg03ZXVy6MuV/wG9uasI72/aXuxiWXoAVDgXCFQ49uHfek7Wa3kAq3b3n+5H7XsNPn1/ZrXNc+svZOO/uWd06RxyaWnvueB5LYcgqHIhoNBF9nYj+RUTziGgWEf2GiC4hIitcBDISqCe3v8UWXC0dKXzyD69j3c6Wol6np9Jd4TB33W78asaqApWmeKzcuh/Hff8FPD5/Q7mLYikikY07Ef0JwB8BdAL4MYCPA/gCgOkALgQwm4jOLHYhDwQyme4JhwdfXYcP3ftqXscyxvDmpqas+xXbJD79rW145Z2d+OkL3ev9Hqj0BLNSKXhnOzdbNa7cXuaSWIpJtmlCf8YYe9Ow/k0ATxBRJYBDCl+sAw9Xc8jT6/B/Ty3L+9qvvLMTn/rjXMz8RgPGDqkN3a/YKbvbu/h039UViaJep6fS1UecOo6I1+4jsrBHsLe1E4vW78XZEw8q2TUjNQcpGIjok0RUr24joksZY52MsZ6vB5cAVsZQ1j2tPNnf/vZU5H65+hwWb9iL9/6sEfvb49mX27t4a9GvjwqHlo60KyBLAWMM9zauxobdrSW7JgDIDDH5doQsufOZB+fj0w/Mw76Y32IhiOsz+CWAV4hokrLu+0UozwFPOZy+HSJ1Q7Zr51q02595C2t2tGBpDJMVALSJhrFfZd8UDufeNROTb32hZNfbvr8DP35uBa55YF7JrtmRSmPh+r0Aes6gyuaOFNo6SyeUy8HyzfsAAOlu+rVyIa5wWAvgGgB/J6IPi3U9eCxw+ejuq8vH9COFQ7ZjcxVc8oOrqcxmfeS4ZqWkuVplMgx/m78BqV5sfunM8966Eyywp4SRQ7c+tRz3z1oDoHSReR2paI3smP97Hif/cHpJylIuZL0qpV8r3lcPMMbYQiI6C8BfiehkAH2ze5iF7n4vXekMEk68R7t0YxOI4CZ9y9b45yp4Wjq5mSpuL0BqDtUhmsOj8zbgW/9cin3tKXzm9PE5laW305VHj1D6OFIlbDCWb/a0yFJpDmf8eAa27+/AujsuCd1nXxaTahyWbmzCkSPqURnSuSkn8tstZYLHuE9hCwAwxnYCuAC8g3xMsQp1YNO9L2ZrUzu272uPte9lv5qNS385Gx0p3ihnGA8zHHfT03hry75gyXIsmtQc4jpa2zul5mAWDlua2gAgtg+jL5GPM1sKlFKaGtQr6ZrDwvV7MKMIEUzbSzDD4rqdLbjsV7Pxg/8sL/i5H5u3HhcUaOxJPp2IfIklHBhjlyi/M4yxbzDGep547QF0tzfV8NNGTPvhSzkd06mYlf67aicA4JE57xrK5hUujkmgpSPlO382pEM6mSDjjHh9PZopivyEg9QcSigclEvpl/3Ab17Fp/9UOv9HIdktgjri+td0Xl21E7946R382jBO5Zv/WIqV2/YXJFihlBFxkWYlIvo3IrrCjLH3FbxEBzjliFZyHdIZhtGD+gEA/vz6elxz+ngcNqzO3U/9mNMZhmQi2mAkzURx7ehy/78v2Ijv/msZnvnfM3DUwf3d7VJ4hPkksjHupqdx9XvG4db3HZ3X8T2V9q40Trvj5ZyPk0K7lGYlNUIpSosoyrUZAxUp7bEsfr6n/5/fv+7+vuHsw33b+lUk0NaVxo79HRgzuAYA18o70xkM6FeR03V6klnppwB+BuAuAJPEb/XPolGO8D7P5+Bfr/di1Pkc4nQ2pQobX3PgwmHJRt77Wrppr3F7PpqDtLk+8Oq6nI81sWr7/h4zQc/GPa1oySPaxnVSltDUoKIKhK0xTaHdobUAEUnrdrZgym3TsTYwir94z3BQDRcAqnnsnJ814rjv5R7ZVkrNIds4h5nirxFAs7I8kzE2szRFPLCQ38sP/rMcdzy7ItYxf1+wMefrqA2b9DmkGfM5ncdrA+KYTzjE/xjiNj5tmtrc0uFfbhdlXvDuHhz//Rfw5qYm1wyWjUI25Pvau3DRz1/Bv97YlPOxmQzLu5e8pakN4256Ggve3e1bX5HIT5PqKoNw85uVvIVV25uLcr1/L97s/t6bZX72OA3nr2esws7mDrz01jbj9mLoJQNrKgEAOxThsLkpP2Ha43wOgh4S1Vx83tiwF+Nuehpz1+7OvrOG/GAWrt+DuWt3xTrm639bnPN1VKdup2JWirI/qymlcxEOnel4PTapGSTFKCnps9C3/23BRuxt7cKlv5yNT/z+dXzqj3PxyT+8jigKKRx2N3eiK82wN8cQ0PauNA791jP45cueRvbR374W+/g5a3h9ePg1vz8o33nHy6ExqNVG/V2IXn0qncF9M1ejvSuNdIbh7hffxpf+usjdni3Zn17fTLy6mr8DaX6VFNMqNqiWaw479ndfu+pJPofBymKCiAZBEa6MsdxbzwOAGSt4xMV/V+3EtPGDs+ztR1aydIahoyv/F5nOsMhGQw3dUwfBqZqDLijUxX1tKby9rRnHjxmYtSxdqXhfjjpCen9HKmAqkcJh7JAavLvLG9U76+0dWc8ttaNc2bavHeu1EcRyJHmuMeNSmDw851387zkTAACv59CBkBqC3qjn6zIoR7oOteRqB6MQZXli0Sbc8ewKNLV14cRDBuHnL73j237lH17Hwu+c51undniyZQgAgE17ecRcWCeqGD6N+iouHB6dtwF7W7swsLYyp+NVTbWUptBs4xwWaMsLld8MwKGFLU7PoBBRNak0Q3OMnkwYrZ0p1FeHO6v2tQU1h7SmOejjGtSP+dy7ZqK5I4W53zoHbV3pQE4m9diOHB3S1ZVcOLR2+u9fCst8spd25PlRnHfXTOxrT+GBC737kykI4go9iWwAk3n29OX84npDGja3987mDgytqwo9X5xAgR37OzC0rrJgjR4L8Vvlm5GWMYalm5owefRA9x03tXUZhc3ulk5kMgyO8vxVAd/SGf97078NuZTPU8pmZpTvd9nmfVi2ORhing3VXJvvIMt8yGZWOpIxNj7kr1cKBkBJA1ER/XimL98WiNmX9STDWE6VVSebmq7mWJH7ZhiQVj4WvdFRhYMUXNN++BLOurMxcH41AiaubVsKVfmBBX0OfDkfLSBf4WAaHCUFqynK576Zq0Mz3Mrn7OTZ0FaI6DC912oanDhnzS5MuW06nntza+j5svUit7VkMPX26fjdK2uM2xlj+Ncbm/Lu9ft6tHme4+mlW/C+X/0XTy3ejISb0I+FNtJ6nVafQXN7Ct+e3Yqz7pzh2+fXM1YFnmNa60TNfof7vnJ9tfvbuzD+5mci9+luRgDVAtFjHNIAXiWiJ4noeiIaV4oClZOWjhTG3fQ0HhI24agcQVua2vDZh+b7bKKA1wCnMizSBnr+3TPx1cfeCN2eLV/Mvjbv3FJQZFg2zSH0dAHUnmDsQXBdMnqG/9fvX1by9jzMbd1Vp9WGTJof9AaNMYY7nl2BS38523gOKeyzhQCH4ZmV/NfV/T+rtu/Hx+6fAwBYsTW8p5ntvWxr5dtfecfs9J+xcju+/Ogb+Pn0d4zbs8EKoDls28edtAvf3eNqZFF+M7VOM8Z8vermjhQ2NjOfyRIA7nx+Ja5/xG8EWbl1v9tBuuvFla4Ji3LUHUzjeXS66xsqtPkuLtmilaYA+IpYvEdM9nM3EZ1PROH67gHKFi2CIMqsJCuR3suUrzGdYehKs9Be8tvbmvHEovBomScXbcKk7z6HpRvNvVhVK5E94XSGBXpEKrmkz1A/9rgNc5sok/wYZBndof9pKRzy0Ry65/BU5VGYWSnbYDIprBMR3Uvd5MQYwwvLtiKjjCt55Z2d+NEzbwHg96U/32894WXJH96/OvRaakNhmr9argozg0lhvdwwmj6MsGilfButoXUikqe5wzUXRc3Frc64eMkvZuOSX3iCPFtDPW+d5x/67aw1uPqPc8X6PXmVHYinReajVZ3900ZMuY3ni1IfR66m0O6QNVqJMbaOMXYfY+z9AN4D4N8AzgXP0vp0kctXUvJJh6v24AGvhyorsW5a0QmzVz4r1OA3Npgrrtp4S9OJ7pAO2FVzCV/NxFNl1+9qxbibnsbiDXtdB3Snojk8v2wrDvvWM1i7s8W9fj4jerurOaR8wiEl1nkr1+5swYRbnvUdc+E9s/CbRi8yyTUrRfgcBmvOxr8v2IjrHl6Av8xd7xMqv521Bpv3tuHIbz+Hzz44311PBAyt984RNkhqb2snvvmPpe5yW1c60KjKKpJwzJ+5NPPF6f1K1HE8vkYrT+FQX83dnruaOyAjelOZ8NFCsu4s3dSE5Vv2+cJDN+xui7zWh+/zR5bNf3ePe22XHJXCOFF/+ZiV1u5swU5RrkKY7/IhpwBrxlgXY+xlxtiNjLFpAK7L56JE9AMiWkJEbxDRC0R0sFhPRPQLIloltp+Yz/nzZXez/yOJ6mnLSqG/LHmIJxyi/Q5htnTZuw4rg+pbkJrDO9uacW/jagBAZcIpmFkpyiEtc+k8rKTrkA1Fa2ca336S94LX7WzJO7XI88u2YtY72SOaovBpDuJ5qQ3afKVXKVmxdT9+8hyf1W5LUxuufYg34lEO6doqf4zHhj28wdqxvyPQ4Elhow4gSxBhzKAadzms8dEjebrSGYOPif+Xvo6OVNo394OMvpKh275GMgay0WrrTGNNzKlh27vS2KbcryzyzuZOtxductDLyD1Zp00+OT0qLS65CEcdk8lI74R116ykPg+1k9Telc75neVCtmlCJxDRA0R0l5hL+lkiaiaixUQ0hTGW+ygizp2MscmMseMB/AfAd8X6iwBMEH/XAbg3z/PnhV5J4giHIH7NIVvEUrhwEBFIIZdRK5w8x69mrHJHYVYmHUMoaw6ag9JwRqmyagPhrfO2y55dRyqT91wXn3t4AX49Y3XkPl997A18JGLMgfosXLOS8gyz+REem+fNl+xQuPkj2DDw51iRoMD9m+pXwiGfgAmrg03agLCuNDOYEb1zAsBN/1iKM34yw+146OM81u3K3sD7xjmI/9c/sgB/eX29sg/fsmF3K37y3ArfM/ncwwtwspI7LOMKhw63/KbnWiVSrkhtz1SXVMGXi5aspjzP1ZtkCmrQ30N3/QRhGtr//G4OTrqteKnKs2kOfwLwKoDNAF4Hn096KICvA/h1vhdljKlGzlp49exyAA8xzhwAA4loZL7XictPnluBcTc9jR2aFI4SDmrdU9NUyPWp2JqD2ezUIT7gsEqeLZ9OVdLxaRdAtC1XR733qEFwci9TZJZa9JaOVE7X376vHR/57Wuuaq0z5bYX8TNlruonFm0KDFpUfRuq5uCOc1A+NN30oj/3/kpYsWxsTSq+focy0qsi4QQGWi3euDdwfNIhnyAzvf4Nu1vxxEJ/vyyVCQpf2eOUms7LYvyO7HXrAibMx9bU1uXZ+pX18nIztXEqsvhf/MtC/KZxNVZs3e9uk/tKU4ss897WLldbNX13Ujh4AsS/fXBtpU+45RvdFuZCaO1M4dqH5rvjJCQmR7yu+XRbOCjPQz2XnHSpWGQb51DHGLsfAIjoesbY38T6F4nozu5cmIhuB/ApAE0AzharRwHYoOy2UazbYjj+Ogiz1vDhw9HY2JhXOZqbm3Hff3mvdN4yf+/0rRUr0dhqDgPc1uK9pDuf9xqphYsWoXldAu0dvFF7dd5CNK/zP2a18s985VUAwKkjE3hti9eYdaR4A/b2qlVoTK+Hzturo1Vhlu7Cps1b8ez03bj3jQ6cNzaJtXvCK6n+/DY3e/uu37gZjY3m0d7vrOMNzJbtwe37m72UCovefAttWXw6ahkeX9mJuWu7cMdjwSwtjY2N2NnciV++vAonVW5Bm6LZqOdo6vDWNzW3utvWb+Ef+KYt29x1K7b6hdv0Gd55GhsbsXmTV/aW5v1obGxES1ewYWhtbXPPuaM1g9/P5tdav3YNunb4BdCNf18SOJ5l0lizdp27vGz5cgxs8puQvv9a0LY++7+voa7Sa9l++JfpaG3vAEDYsX07GhsbkUnxe5gxazaG9HOwcq0/AOP1efOxY4BfQLR2MXzhpVZcMr4CHz6yEi0tXu983/79xu/u5cZGVDiEPSJF++tz52Gbdt4XXp6JmgrCEuW5L3uLp5vZtmMnli3b6z9pmu/331dfw9B+Dt7c6X9fo/ulsWSn9/288PIs9K+irBqEXv6mvXuN9/Tq5hReXN6B/Xt24vrjvCCB1XuDHafGmbNQpWii+5rDzV3Z2q1xNz2Nn5zpjeZ+Z/VaNDr+jkFzc3Pe7V8U2YSD2proIQ2R4pCIpgMYYdh0C2PsX4yxWwDcQkQ3A/gigP/LVlgVIbTuB4ApU6awhoaGXA53aWxsRHVFO1o70+ioHADAa+QOmzABDaeOMx63Zkcz8Eqw4TruuONx8qFD4DQ+D3SlcNiRR6Nhsl/5ae1MAS88DwA4YcpUYNYsNBw/Aa9t8XIx8UhvhnHjD0VDgz/LIwAsTr0DvPN26H3V1fRDoq4WP17UgnW70jhq/MEYc0glsMZsntGf31tb9gGzXwEADBk2HA0NxxuPWz17LbBiOSpq6oA9/irCklUAeAN08NjxSG5cB3SE20jVMsxoehNY+y6OnjgBf39neXC/5552f89btxuY/lrgHKt3NAMz+DuqrO7nbvvJ4lcA7MOgIUPR0DAFAND+5lbgDS/c8ZT3nA688IJ7zl0LNgJLeZqT+vp6NDSczk1mL03HxBH1bu+4qroap51xJpIO4YJ7vBz+kyYewVM2zI9OaV1VWYFRo8cAa3inZMKRE9EwZYxvn7vfnA00+aPYTpw6jQ+Ym87LfP+SDlx9dBWATow+eCQaGiajavZ07O/qwPFTpuGwYXV4YO1cYIvX65983AmYMs6fEeDtbfuBl2ZhRXMlGhoaULOgEWjhPfSa2jo0NJyB5AvP+LSd004/AzWVSQxYOhvY34TjTzzJHYmfeOEZpDMMJ518Kob3r0bLki3AG3xs7aGHTwCWL8OgwYNx9NFj3PUA0L+uH/Z0tGLqtJP5gM2V233P8qIpE7DkOe/7eWhNFf7++fdwDeV5f6CBilqXAGDQoEFoaDglsN/uhRuBJYsxcsQI37dQu243MMdvznzPaaf7BrAmX3sJaDOnzjjrrLMw8+0dOHPCMH+gg1KmE6dMA2bxenzwmEPQ0DARj85dD4AHJNTU1ga+30KQzaw0UTiGlyq/5fKRUQcyxs5ljB1j+PuXtuufAXxQ/N4EQP0SRot1RUWqrHpmSV293dnc4drQwywkcrU89pV3duDwbz3ji6rwO5X4b92RKe3hYaYY3WSkk3QIq7btxzoR8z2otjJvh3ScCAmTg1BNZ9DSkcrJDtzqTlEaNHXoz2RnyGQw6nM2hrIq96X7mPUxJh2+RIcZ8Z/vc81p43HlKYegf3USjAETbnkWN/59ie/+ieLl70k4fl+R6f0fdlBdYF0qzUKjlaQ/RZqX5L3tbe3ypWhRn1cqncHe1k7XeS+jpvxmJSbK7H94svy6ExnggRKAOnAzaDKJMivJc+smtBMOGegrx/x396ArnSnYXBfyNLrZyWQy0j/NKIf0M0u34uo/zcODr60L3cdnHhXv6KYnvEi1YqXYyiYcJgG4DMClym+5fFS+FyWiCcri5QCkyH8KwKdE1NIpAJoYYwGTUqGpEjOX6Z5/tZJu29eOKbdNdyNWwpyrus/h0XkbkMowLFq/x93H39DwFy9D+nTUurdpbxvumf42GGPoylLpHYfcLKgAb2RyckgrNTwqjFQ2EK2GkF3VGf/rGauxszl+VEirHKVumL9aF1ZhIcjq/fpCWeUI6bS6Xe8I+MvamVJNfhmxjv+vTDq47f3H4oKjR7g26b8t2OjLtppKx3v+RNroXXHM4/M24GP3hzvcTdFK7Sm/z0E2nrKxae9KY9RAz2TR1NaFVdubwRjDnc+vxPHffxEfEuGf/aVwUB3SUvhowiGjCYcMYxh309O4Z/rb7hScMrWK/x15Pgf9UclvNMznUFeVDMwTsi8kDYeKHmYa5nOQ5dTHNZh8Dro/MKoMm0V92bgnPAxXT5+h1/dipVuKNCsxxoLTiRWGO4joSHDT1LsArhfrnwFwMYBVAFoBfLpI1/dRJdJk6KkWfL4B4Uh7Y8NeABHCAeaejRqvrg6Hl7/rFM2hKum4DVCaMXSlM3hQzGNwz/R3cOUpY7PGTicd8vU4UhEpCUykQ5xgYexpzT8c0ETUgDM9oaE+1kSivoJORROTQksVMvo97mrxOgp8JK76zvzjOWSvVi3qUSP7+4INutKZWJpbVdLxD24Tx9z4D88/0dGVweEH1fnSZHelMwHNQQbhSGe7bKw/dN9rmHPzOUhnGIbVV7khoL9pXI2lm5pw00UT8dwyf7oJ03gLt8EM0xzEA5HC+J7p72BYPR87K9+v+o5krzidYQFB52oOafP3VZV0kEw4ALxn3tTWhUE10Unubnv6Ld9y2Ahp2QnStxqjlbSyRX2rsr0g8AG1yQRh4oj+vn1ULbYrnQmkR+9Gfs9IsmVl3Y/omeD6h22LgjH2wZD1DMAN+ZyzO1SFzEzmC38UFbyuKonH528wRLMksa89FdAcJGokiBr9I4VAbYhwyGQYHnx1HW57+i2MHMAdYRnGssZOO+QXDi8u34Z+OSQSVBuoKM1BfqQdqQyG1lXBoe7P+fv8sq1uZI3p49Mjo2R+K9005EvvkPGOla9GvUfdlNGsdBR+/NxK3DeT+2ouPHqEO9JWCinZG1YblpEDqn2RLV0xNIcjh9ejuSOFdIahIkH8GINEaetKB95lytCgSke9NCupZpd/LNyINGPuCGXAyxDw9rb9OGSwP3Oua1YyzAmiaw7yWcoAMGkCGVxb6X5rRrOSojno9y07cHJ//VlWJBx3PIekqa0LdZpGftzoAVisZB3QJ48K0xzk5eT2Vdub0bhyO0YO6BfYNxezknpembZl3R2X+PbxaQ4phtWacCjWNLHZ0mfUCwHwcwA3gUcOjQbwTQD3FKVEZUCqrABw+uFDser2iwCY0wFnGMONf1/iG50KAN+8aCIA/rIzBrVY/XDV3EKyAfdpDsqHn2bM1Whk6CFjvNHUPwaVZIJ8PdX1u1uxctv+0P11ZA+tvioZOVZDFRynHDo474ylasN1u9KbM6ntenm8EeLhc1bIYqp+AJ9ZSbuOeg0pGADg4IH9PLOSqznw96VGw3akMq7QALggyuZzSYqxEKkMc23zJvt7W2ca1VpSyK50JrBva8pv3lFfTb+KBNIZhuqKBD516li+vxC6+9q6AmGusvOiXmFPaxc6U5msPgfpbzMLB+W4tPTlBH0F8hl7PgffZlQmHTfrraSprSvwXj931mGIQ3NHCq+t9oJT5PWkWenPr7+L255+C3dPDwaFqB0axlikz04/r4l2zaykp58pllkp7gjp9zHGfsMY288Y28cYuxfcV9ArUDWHigS5lfpnL76NP8xeC8Cza4dlS5WNIkOwBwf4Gy3V5yB7BarjtVKxVavHpZWeVTrDAk5slXwnkJHICj58QDV2RfgKVOEw4aD62KmhJ46ox/NfOdNdVhvOMYO93pipcdQz4arpy1VfialXqtprpebQ1pnGa2v8obhh77mqwnHNRbrmoBod2rvSPkGZimFWSiYcLhzSGbeDYNI22rrSqK5I4NozxmOgmIKSaxn+/VqFfKtwpEPaq1e1VVw4JIhwvWgwWxVH9damdnxUiZKSxVCLs7ulE1/488JAXZN1Vm/wBtdUolI08m1d4T6HjlQwFYg3zsE8CI6blfj1powdBMAsHOJ+F6f+8CV8/HdzXIHpmn+0w00z4KnvQR8XEdjXVR3C91HNSm2dqeBYmjILhxYi+gQRJYjIIaJPAIg3Xv4AI5lwfA3cD/7DwyjD7NoSaddlzNyg+eZHUGzRcpCcqi1UKb3CdIa5X6SqUnelGWoqEqFqcFhyuP4hjm8dqQoP71+FXS0dob1eNbVGTWUi9seXTJBP5ee3yWf/+u8qr6E2TcijT+qiNviq4PWNLBWr5XscWFPhCoebnliCf2pJEMO0Je4T4MJZmgdNPoeOVMZXls4YZqXKBCGdgU9zyDCGP/13rbsPYwztwqx0yyVH4aFrpgGQwkczK3XJBi1oVqquSLhzI1RqZtVdLZ3Y2dyB4f2rXB9BWNmnv7Ut0GNPhQiHQbUVikM6OMhTvo9dzZ14UpvCVfc56N9YTWXSDQA4ZAhPP7KvPRXotUclTQT4s9q4pxX7xfuXHQAvWkmmIglvkdXO4ek/nhF5PbfDF2F60jPP6oKzLGYlhf8B8BEA28Tfh8W6XoFagSpD5vPVe6s60sSjJ7+TpBmfe/jNTU2+3rZMVleRIEwa2R/fv/xon5krzbwkZPK0jPHGIJlwUJ00+xHCGumh9fGS6aZc4VCN9q5MaE9avZd+lYmA3R8A/vOl0wPrkk7QRpzOsEDOINOz1NOcqIJ78q0vuAJXbXgeWt6JXc0drpYxpLbSFYBLDfM3hI1sl++mM5UJaA6OTzikfc8sjkOaj6Lm9aciyU/2l9fX43v/9sZ5pDM8TbVMJ+9NIBTUWFuEWUmuTfg0Gb5/gigwh/VakQdrWH0VXv7aWe51+bl0Wz8F6prs3evr66oqUCWu9cic9Xh5xTajmW9XS6c7nadEj1bSZVV1heNqakNE8sN9bV0Bn1UiRrp1NRNySv3o4L3jqMzC33piKR56bZ0vnUcY8n50M56KKhxaOtIl0xxidSMZY+vQi8xIOmoDF2bHz5axNeGalcySPJMB/jZ/I278xxJ8bKqnrqtROc9++QwA8KVGYMzrtblhfCKUNZkgVFU4vsqjlwfg9mW5T/+I2eVUXLOSSBm9q7nTaMZSn11NZcKYsXTEgOpAnH/SIVQlgk5Vne/+a1lg3SYt7E9/N3taeVn1063d2eIKlhEDqrFmB1d+TZ3iMOEgBUFHKh2MVvKZlbhtePSgfti4pw2pWD4Hx9UKZWO4TpubIM0Y2jo9h7Q3gVAwWqlFPBa5Xu00dwofRSJBoQEZ1RUJ1FdXoL46GdooD+hXYRAO/L++/h8LN7q/39qyD9c8MN+3PSoljNSm05oWLSEiEa3EtYiqpGM2K2XTHAAs2aQKB7/m4IRoDhcfOwLPLOURXq+t2YVFG/aEhqeryBQqeyOEQ3unF+6uBlRIipWoNVvivW9r80jr299LRJcWvlilRX3RyRDNYV9bKvJlyx7cVx59A9v3BUdDphnDm5t5pVPtlLLRVlVz9WM19ZwzUnNwKJbmoDovwxoCHflRjRDCYWeLOQKpQxMOpo/PtD7hkNs7lqgfvEnzkfehxoSn0pnA4EXZxugNiOOQex/D+1djS1M7Vm7dbzSZ6PNfS9TZ3ALRSgbN4YKjR2BoXZXPrKT6VHzndsidzS9Mg2XM8zkAXn01jXNo6fQ3pKog70xlXJ9D2LXkfSUcL2mg/qj6G4SDbFCzNcQ6UVE9+iA403chr5ZwCPXVPJBCD1HOFjBBxNPQS/TQWXm0Hk79oysm474rT3KX27syuCfGJEqNK3mIfJOS/E/vRDwj0vf3r64wDiYtl89hKYB/E9FLRHQnEd1IRN8loofFKOnLwBPyHdD4NYcQ4dDehdFKKmUdWema2rrwEyXXkiSTYW5FVV/t4/N5KinVbKvagLk5SjuXMD0kHcfnn1DR7cum31Gs2sEF2PD+3AylpzOXqM+uuiJhjLqoTgbXJxPBRknVHExCTKay3rjX+3hfWbUTe1u7cO6kg9x1cp4CvdHnKY47UVuZcDWoC+6ZZWxowjQHKcRTaeb6W2QvX73Dts40WjpTqK1MoDJB3CcgHtX33ne08dwVCQeZjIhWChHi6YzwOVT6NQdTVlaZckqu9o0Yl5qDQ3AcMjaa8ltIEIUmoexfXREayqrWwbD7UYkaT+OaldJmIQV431DC4XWr0xD1FDUXB8DfoarB6BqT53Pwdx6SBvOaPiNdFNKsVF0RzKYsO5N1VUnj/DDlCmX9F2PsNPBBassAJMBzLD0CYBpj7P8xxrqXaL8HoL7oyhCz0v72lG80qY6a8tk0LiCdYehMBXs9Mm2yX3PwGvBMJjjxCROmh4oIk0AyRDjE0RwWrt/jzgshJ7gPs4n6zUpJ48fnOBSIxkg4TuBjUu24JiE2qLYSNZUJn+bwytt8CszzjhrurnOjO7QH19GVwa7mDgyu8w+MMjU0YQ7ppNsYZ3wjpAH4AhnkmJeaKv5M/rZgozuXs+7AVc/No5VYqHmzM5VBV5opZiVPWIVZZWRPU/Wt8RTqnpnE1CmS6xwnmG5cYjIruQ5pzbRp+u07LkJzcO8z4+/Jq0izXlI42fmzylVz8IeAL964F88u3eJWJfmKdc0hmTAL2LhI02j/6opQQSzNSoGU4GWOVqpmjD3AGPsRY+wextjzjLHoGK0DiDhmpbbOtBu5YUL94GUj99nTx+MqEUMuRzoD5jTc6gcWiFbSyDDeu0kmnFBNIExzqIghHNTGt0akrwiL1e7UopXifh9JhwJhr3taPAFkahyrkg4Oqq/ylW9ncwdqKxO+RGfSVCcf3Q8u5z319q40drV0YkhtlS89u9GsFCIcfGalVHi0kmTMoBp3XItMzhc2d0RFwhFmJRYaUCCFljSxyQbJZFaSqGalQ4fVuvvzOsSPN/Xs5b2qmoNed+uqk4H7dtNnKOvVUO1jRw0wljNKc5CfpRvdYxIOJPflTvbOVCYgcLJpDoD/Hr/86Bv4/J8X+sJzv/rYG5i7brdvbFKF4xjPTQSs+MGFOHfS8MA2FW/ej6DmIOnfrwKMBcOsyz3O4TdENJeIvkBE5jd7ANORxayUyfCBLINqwp256sf8uphXYNLI/vjYtEPcc7iV39CK+IRDDLNSV5pxp25IYy/PR+TXhgjA999Tjc+cPt63v/pBqD07+VGHjZLWo5XihrLK/Z7/ypm4SQwgVHNbmXrXVUkHw+qrfAJzV0sHaquSvl6bJxyYKBf/iNtT3Kw0tK4SO/aZhcPsb54NwD+969ghnjnRMyspmkMi6JCWnHnEUDfdhiTMdFmRIKQZEwMczfvI0eGu5pD0fA5hvfuX3tqO1Tua0ZnO4PTDh6IiQehMcTOX1BykcDCNt0k4PMQ2nQkO6CIEHaJdaYa5a3dj9Q4v2r2fct6aqhDNIcI8IjsSnkPasI/4LzUHKQBV4jikTR0yWbY9LZ3u3O9DFA3Ucch4bsZ45yzkdbrIaxKFp9uQPk9dqy1rKCtj7AwAnwDPmLqAiP5CROcVpUQlhjGWNVpJmlSiBp2Zjks4nh1S1RxM7ae6LmhW8r/8TEaGslKE5uCIcgXNN4f0T+A4kULZPadyCbUXKT9kk3CYsXI7Zq/a6S73M/gcwmSFbMyPHFHvCt1dSoiqqXddlUxgQD+/SWjnfh6ZpGpE7VpqhlrRMLV3ZbCzuQODayuxfb/nxFa/Rdkgqik6ThKDqgC/jb8jxR3Hsseo3+thw2pRX13h+kD0e9dJilDWVITmsG4nt2PLKLIKxzO3hGXwXbOzBef8bCY6RXllrzrNmNtoyfv29YZliK7DzR4fu/81X0LCcUNq3EF7KgvX78FHfvuaL0RYFTph9x+tOYjvKBOufctOVyLhcLNSOhNwcmfrvBCZBY8cDKeOv5Am1zjnznZd1U8T1tjLd9OsjfMpt1kJjLF3AHwbPHXGWQB+QUQriOgDxSlaadBnwDT12GQStqh0FaaXn3DImxc3w9zkb3rPJKGZWPyaAwK2c59DOkxzEKercMjYC9eLK8v0yJx3cdUf57rrayPMSn8Uo8clJrOSq8GIZfkMVfOdFGTq+AXT86xMOujfzy+gt+5rR21Vwm0kAU9zkO2H7LW2dKTE4K5qXPWece7+rYogkO9fNSupjYCrOWS45qAKUr3jKI/TGyi1jh06tNZbL6KVUkIrNLF0014AwLGjuQLv+kBSwfQZOrK8suGU0UqA1yHwCQfFIf3i8m2Yt26P73zVFQk+XkK7rmlUcE2Fd96whjIqh5d8HlGzxfk0hwTPT5b7CGmzf0XWTfWyQzXfVdS5o9JjAF5dTRCF+l6k6TSoOUSeOm9iCQcimkxEdwN4C8B7AVzGGJskft9dnKKVBv3BmoSD1BzCVP2wbUlFc8gw5qZ9lh+BjATSK5Xa4GeYySENzyGdRXNIJhxfJJTrWNNMIBmh2dyj5YqR5gvTh3vUSH/eRZNZSV+WFdwftsr/q7ZUU+NYlXQC4zSa2rpQW5n0Ce6VwrYvryH9Jpv2tiHDgIP6V+PTp43Hjz5wbOC6sresfoBqZs+kT3NIa8LBX+YhdeaMoKpWpNYbOc5BCn4TSzY2YUhtpRti7PocMua0LSqpDENFwuENZ5ccqMavIzWHWp9wkBqRuWGTIa5dmQzOPnIYHr2OT5Lz17nrA/uqzyns3qJyEKmdLCDErKT4HKRDWjcrxXEaRwkHlcG1unAIP2dWzUGGylL4eA9pVgqk7M7y3vMlrubwSwALARzHGLuBMbYQABhjm8G1iQMWffyYSTuQEUVRwsH08h2H3J50OuP1IKX6/OGT+GA4vcKq1+FJ/DSzkrBLq5qDHgEiT2EyKwHmDKaf/MPrgXkMHIe4jdrw4errKrXUI0DQxisruGqHlw1Uu2Gd7/xJx5g6urYq6dNEfjtrDZ5ZusVtQKTzVsavDxeBBSYfkufg9Z75YcKJe8aEoUp0ENccqiI0hyG15gAG9f2qYz3kPNNd6YxvJO+phw5xf2/a24bRg2vc50zE348aKhtFZZKblaR25ZqVDJqDFBhhTlxpAkmnGYbWVbnHmtoqtQ4mHMLJ44PDp/QIINPx0uRiMivJK8ixG13pDJ5d6k89njWUlYJZVQG/cDDl0uJlDG8fsokkNeVImAYoTXP6ALxyO6QvAfAXGaEk8ivVAABj7OHiFK00bGnxP9kozSFssBBgFg5JxayUUXwOnWIAmzxGb0DVU5kmPpHhjskEuf6JfpW6cJA+BzI6SvVGPMMY5qzZbbw3GTOu09qRdrUfeU7TYDcVKRzUzLTyGNU2b/qGq5KOO+mMSm1VMvBuViiD2xwiVDjAuyKdwUGi1206l/r+ayoT+MNVU3DeUcOx+ocX48FPT/PMGxkWyLyqP2e9Z+leQ2lEfIJCCARZP9z7VqLXTCanpJg9TtUcTjxkIH5/fk0gwq4qyTsUUhDLxlLeR22IWcmEbMjkaP0o04la5qRDeOxzp/JpUxX0sQP+a/H/avJJU3kAMYYm6WDT3rbAvBRZQ1nBvwW9PqnzlQwWmqT+TUQ5u7MlpFSjocIGA6oDHlXK7XOYDkB9kzVi3QHNko17cf+SDgytq3Tj5E2OUFdzSEa8fMM6RzUrZTzHd1eKhQ48ksdJTGalDOMJ6SoSDqorHBG6pzcYfLki4QQc2kCw8Y1KCidVdB0+yCuJEw8Z6K4LMyvJb6O+ijfIakPgzlCmmHdM3xJ3SAcb9LqqROC91VclXaHqEKEyoWgOQqCZnPlqEEF1RQLnTBrOhZ4cLKaNc4jSHC47zj93uMRnVlIEhZcnKeMzvfhnlMsE3l3SjT7y3qFDJDon/n0rEn7NQdYT2euvUyKJ3JxRYQ50YVaSZrCwEF3AXy8cpW6qtHWmQwfLyalRTWYlqVnpZiWTJpLN9s8d0sFBiKrmIP1eujCLUByyIgU7gzk3G+BFHQZ8ltktZXkRL0UnH+fg5nxgjDVLzeFAhkeFAPd+8iT8eQ63k1YY3nAcn4OpzqnagRqt1JXOIEHkmg70qqBWYFNFYYyr8gmHcNzogTjl0P1Yu8OfJFdeN5kgo5of1Bz4gCbTYDcZFqjT2plGTVUCf73uFFcTkKeVk9V4Dmn+X35Yqubgzm0ckcxMlkPNKltbmUBLZ1r4HPzvht+3HIzFnb0tXRkQeY7iME2wIkGiwTNpg/yYtBD2aiMidz94QDVevfmc0PsIMyvJxrUz5dccfPNCZFjg3VUmHJ5bSXnRsgE2aaWVSccdKCjrmkz9rWbKdTWHkGrvONx52pXO+ISqCVVwOEodUdnc1I6JI+rd8SCSv1x7MiaPHgjAn18MAP75hfdgwvB6ANoguIQTiBLTy2GCwIMCKpMOoGSMUYWDO/YnpfszwtuHbHm15GY5hklSXeG430qFqzl45/rQSaNx6TB/oEChiCvrWojoRLlARCcBOOAHwZ14yCD86PR+OGnsYK93ZKg78YRD8EA1WomPc/AGIyWjNAfVrMSCNlyuOXCH9PtPGIU/f/aUQO9OfqhhDWDQ58DcmeYkdwiHbUWIWamlIyWSnHk9enlddwKcgM+B7+f3LwjNQREYpm+pSvE5VCYc91w1VclAQ9PS4SUoc4ggO8RDaqvc9xjWS5XbjUEG4jqvrd6F9lTaF3YsG6dsJoQKg0NarSudKb/PQX2H6Qwzag561JDcRS+L9FHIhlM+e+norzX4HMLMJarmUJFlhLBaDvluTY3pkSPqA+vGDKrxmfMAzwxzzKgBnp/E1RwcVCQds+/DcC//uuE033I6EzQrqX00eb1AOnDlkINiZj/WkSZjU3llfVU7arnmr8qFuMLhKwD+RkSvENFsAI8B+GLRSlVC9FmyTLHicXwOpleUIK839Z1/LXPn6u1IZ4TJyXw+9UMyjnMQseXqx6U3wqrmoPK1844w7p9hCAiHD5w4GgCvlB0hmkOt5uuQ55XmFvlRu2Yl6ZBWhI0+8X0YlUnHjc46+dDBbq+bm5X8z7K5I634HABptld9JOGag+evCdv2+9lr8d9Vu3znkPeof69faPDPPqaWVbXry9VdaeZO0AP4xwh0pTPBPFUOD01VNQe5i95jdYh8moN89lIoqCGn8vlmc0in0gwJx8lqspHIEplG65uEgz8k3JwlFQiGsoaVWUcfsc0YizQhDxBaVr+KBM4/ajgmijKr3/PMb5ztOybuJFhg/kFtannVYAhJd0xZ2YibsnseEU0EcKRYtZIxFp3D+gDDcxwHt+WrOSQT5lGTXekM6pRRvaYPWJLNIS0Js/XLCBgAuP2KY3DosDqsR7ABS4swRxXZOJoc0pkMw772LndiFb3sVSH2atnbN/ocsoRdVCUTOHbUAPzvORPwqVPH4vJf/RcANxNVaNfhmoMMDyRhr2Xu4DEgPJWIaSyGRO8dq85iL4LIf8yNF07E5NEDcP0jCwPn8CKCvGfXmc74GppzjxqOBe/uwYqt+5FKs0AjXOFqDt46qcUETZa8TrS7oax8P3ckvDK/eRyHNE8UyKesVevj8P5V2KaMQlcLIt+L/s4APo+2TkLRsvXcSuopVJ9DtswBvvtQ1jEwblaK+Nbfc9gQHHPwAHx06hhf0IH6nGSE3JCQoIQw9EGFCUNdyTZ/fKHIRe5MBTAZwIkAPk5EnypOkcoDucIh+OD3ikiFigThyRtOw8XHjggebzhnwnGMkp0xv+NTR11tKg9THNL6MZNHD8CL/+9Mt6JWOJ5DWm2U9AaGMf5RqDZ9Uhr6ptYudKUz+O3M1Rh309O45sF5eHdXa0BzkEWSPXy9Me1vilYyaA4mk0Blkjvfv3reERhaV4VbLpmE/z1nAi6Z7Dl+h9RW4uAB1SK1sXevsg1X1f1smkNUplLTOeTeph60us7fG/R8A2oHQ21oqxIOvn4+75elMpmA8BkxoBpvbm4yxsfrz5GIv083lJX8wkFNGyLvP8ohzVORi/oszlVflcQVJ4z2l0OVDkzee/D5DzI0pjIYgIeZesKByN8pkwJROqQlFx0zwneuKJ5ftg1LNzW505maqEg4+HzDYYFoNPVbJyL8+n9OxJPCZJXN5yCRubVM5fV8DkUKT9KIJRyI6GEAPwVwOriQmApgShHLVXLkOzC9RFdzSDo4fsxAfFWYZvzHByudalYybfPmnQ4/V9owzuE/S7agQ5vYXR7TcORBmDC83mhWUs1QenEzjN+7HhIL8Eo5d91uXPfQfNw/i2cWlXno9fDN7JpD0rcd8D4APQeRjt4bvPjYkfjqeUegKplwHamfOPkQ1FdX8OkUlRz8B9eFx/HrVLpmpXCfg1smJeLJNSsZzul3yKpCwHtO6qlVwaR2JLoMmsOHThqDNTtasGj93kBZdEgzK8n3I5+LOlrcnV40THNwyB3YWZFw3Hm6B9ZWBJzYqtyS78XkHDZlbJXXTzrk0xz05+CZ0vzvrtJQ17IRlWI8LFOC7kO5ZPJIjBkcjNs5cng9Tjk0OM4D4EI03Kzk156KTdxopSkAjmJxxd8BiDeS2b+eR/DwD8Y1ARg+FtP3o9pKTdtCNQdlfcqQsluOQDWp5bJR8Xwp3qCasAYK4B9bhjHjxyk/tBkrd2DskBpfDqR1u/xRUnrcvOtzENv7VSZxy8WTcPZEb/4FWdY40Uph1FQm8fZtF6EiQZi9aidaO9M+u/TxByUwY0PKF76qn082kK7mYGgE9Aag0qe9hTc86jb1tcmGT/aOvev4tQjHMR8PAEcM52GeO/YHJ2TSa49DPGOp65CWmoO4d9MkR+o5PnjiaHdGt6RDru8o4RAOHlCNz515KD46dQye1ObkVs8h34tJONcYOieOUqfVUNZgJJan/avvNu47UqmKMCuFRSVF2f9VDaciSbjmtPHGcUWZjD9aKeE7rrSaQ1zh8CaAEQC2FLEsZUVGwVRrk+cM6FeBLU08MMszN8SzxhGF91QSDoWG1amHpDIZo4kF8NvE9d6YPr5AXaevl8dnmDn2X43K0GfDC6jWrt2dfMtuGYhw7ZmHBtYBWrSSVoahdVWYNs7c25KoA7n2t3szZhEBk4cmcN+VJ+KMCcMC+wPAnR+ajIuO5eYp+QxNYc16r9Hnc4gom//ZKw29EnKqrld9DtzZmzYeD3jvTE0WKNE7OySu5UYMifs5QSRivPiYkXh6if8zV8+hjip3FOEgU7DffPGkQBn1ZxapORiEQ1J5Rmq0kl6HXc0BfoGgvqO4cy5EOaTDNIdcIofChBTTopXUToG8J3V7MbvrcYXDUADLiWgulOhfxtj7ilKqMnDD2YejvjqJD57ot5XWViVdB5Cba8YgG0x22VSahVaYhBKtpL9gteK0hUxXCfg/LnkK3cldlUz4PmDTNWQZMowZhYOae0gOYgOAmy+aiI8q82EDwQE5soyesAjeRxzN4Y4PHGu0R5uoq0pia1O753MQDdeFx/gHpakN9kH9q+NpDlE+B2lWCjExmpBFUFOt6NdOOv7ABv0ZyqliWzuCgwgDdcvRhJP4PWZwDdb+6GIQEW74i1bIkAYoQeSaAvVGVz7bgTUVmHXj2fj644u900X4HKoMtv6EUXMImpW88/s1B/WcCYfw0tfOwoJ39+DGvy8x3xiiHdKhkzXJfGaGSh4IOgk5fUaLVlLPVWqfQ1zhcGsxC9ETqK5I4LozDwusrzcOCjKYlQzn7BQhq5LjxwzEGxv2uueIM85hX3vKOMIZ8Pds3ayOYp3MkXTFCaPw4GvrAGg+B+1c0qxkurc2pUeqDpL65Klj3QFBXtmDmoJK2HgQIDqUNRHSWzNRU5kUZqVgREsYVQYzRL7RSiaiNEhA+KfILxDUfdR6FNQceBn2K8kCPV9Q0Kyk3laYRqOi1j/1bKpZSX9W8rQD+1Wgf3WF7zh5PrNwCK7z0mI4frOSY65bDH6tUE3WSEQ4bFgdhtVXRQuHCBNmqMYvDonj1wgVbGC+lNzZfA7F1Bxi2UcYYzMBrANQIX7PA0/E1+tR0zWoE67rmN61bqL65oUT3d+qs1pv/NWPtKm1yxheq5fDDQ8UFejzDYfhh1ccyyOrpOBQfQ7aPWQYt3eaerh+zcETBrpgUM87UYQkfl6L8Q/LQQX4E6/pPa1cVHY+OjtjjIUPQ22UZIMf5dNx9zVoDibCI9O8XnFYRJPun9JPJZ3ib23Zl7UsRKT5P7I/m7D65zheQkb9uThu3eaor1Oa1KPGkaiofjTZMKZNZiX5QwtFHdAvWE+rI6KRAERGK4VpDqS8y7BtkrDnnmHAhj3e3NO+UNYYszgWkrjRStcC+DuA34pVowA8WaQy9SgGKsIhKu5bj9q5+aKJOPpg/+AaX6qFSM3BW9+ZzvhyDqmoH1cq7e+NDe9fjf85+RB/eGREAyM1B1Od9Zm2srQlsic5bmgt1t1xiWvKkYeZblk2JFEpm+NGmch95f0A0Y22xGSGMPUQA+nVK4IjpE2EhYPKd83NPd56tZHUBYfesOidEJVAKCsQWifCz2GWDklH9V34yxD8RlTtI1xzSDiE68/ydyjks0s65I4BYAazkivgyf8+Bxiy70bNzQKYOwaSMM2hWlzTFM2oE2USe3dXq3G/qHFWxSDu1W4AcBqAfYA78c9BkUf0EgbWxNMc9FWXHz8qsE9FwusBJp2IMFftrew15DsC/B+knBN5WF1w2L78GKNMCBkx2M5UaVsVc88TC3kUyhkThhrLtFtMjKRPhCIxNZJhDZQ6lWncKBPA62G6uZViaQ5eIy81B5NZSX9u/QyhrMYyZfE56GHP6m+9ruj3E2Uf18fJOOT3bcQRulGag1pG/TqAJ5x8mkOEzwGAO22sTjKhhrIGy37bFcfggrFJnHH4UN+5Tckas41Yjuq8hYayJhysu+MSfPaMQwPbggNdzddlzMseDJjHOZSKuFfrYIy58YtElESom6p38KlTx2LM4H4YqEz0UqFFAvnQVpl6FxUJx20keOUzO6T1iru3NTjRCOD/IOUI5hFaCgz1/GoDpZcuIxzSJo1ZT/532LBaPPDpacYy7RK+Dj1VtLygqaEeNbCfb0Y0yXcuPco9Ty6agwzflcWOc6TPrJQMNyvpqA58ubvpqGw+B0fTDgI+B+XwYJSOt2LUQH8abNMIaVOGVJXDD6rzmQ99/gJlQS+j7zrasukc2XrvOknHSwCZZsEcUwfVV+Pjk6qQTDi+svWriOta9dCfk16O7sBYuCaZYQwbdre6z1/vMOTQR+o2ce9yJhF9C0A/MXf03wD8u3jFKj/fv/wYvHLje305/6UD2OyQ9q8zSfmKBHlOq4jKpzegYZqD6RpqeojAeSN6nxnGxAcXLNPnzvL3hPiE6eayy+yV+vy6ElMPOplw8P4T/JoW0xr2XDpNSYendcjkojkophnZ4Js0B51+laqTP1w6hD0v2bAntQ9fvXZSy1sUdT/6HBWmjofPt2E414v/70ws/r/zw09iKIfe0MvblVqr2nP2Qllza2STSrQSY8HstCphWlhcEg6F9tSzZXY1Edfn0JHKYHdLJ0aJuS784c1BbSYsWKUQxH07NwHYAWApgM8BeAYH+AxwcVF9Dt5E8gbhoK0y9YqSjqOM9vTy35t6dypyPgkdU6WPm8vFPAiOr3/pa2dh9jfPdrfdfNEkfOOCI93lKPVWTooS0Bzc65qPG2iwCwPec83PrBT/2KqEYlaSmkOMRkB1bObjkJaP0qFozSHKIa0i05O4kTtaw04UPvbF28cfHaUqjudOCg5e5L81n4N0SEuzkrLNC54Ir0d3feQ4nHXEMN9I4mTCccPKw4In3H19ZQvdLRQ5h4eJ7pp3iMLfoQz+kJ08dbcobaYYxE28lwHwO/HXpzA1WiZbZDAZmiHqQhnp6jjxNQfT/LX8GoZyGM5p6lvo3xVj0skHHDasLrD/NaeNx53Prwy9rkR+vGGzoIWp0ya7MODPlxMXR3dIx/iWjZpDDPNBtc/nEF7GMAGlRivl63NQ0eurKTVLWFRUGLJ3+ttPnoT3HO75mvRR3Pp1fOdQCiLrRpTZ7gMnjnazAksqEuSOHjaZlVTiPq8wHIpKk5H7+YJCOvocavZgr0zSFN2zciutJaI1+l+xC9cTUCeXl8QZ5xCWtM1zSDuhlUytN1H1UG28Tj98aMDeLHFHCodcA/CilcI+pH6VCVxz2nj3PsK49xMn4vyjhgf2kWcNO78uHGSDlI/mIHPw5KI5+EbURgQe6PiEg/ZfL5MJWbaKhOMfBKcKCm0KzqiGRT5Hd4/AAEu9xx96Khdv/gVNAOTikBbrzztqOL73vqMBhGfFDYNHK3mD4OKalfIRDlF50bqrOcjEm1GMMJiHqSdqDvAn2asG8GEA0bkMegmTRvbH188/wj9XcAyzkjGFt0M+h7T7oiNGSH/zwon40bMrjGVTe2uPfPbkyPvQy2iazyGTiW54oqJ4JBcdO9JNQ+G/tjTJmY8baBDCgOpzyEFzIB5imXZ9DjGOUXaSDX7YdI0q1b5BcOH7ZdMSKxOOlpXVn/LB8b278OvIAV9krlpwiDSzUvaGTj4GWaar3zMOdVVJv48kYFbyn0N2UGRiRH5Mbg1dMuFgxdb9WLR+T9YGNo7PYVBNhW9yIxUiL2Dk1EOHYOu+dqzd2SLKkXsDHfQ5RO8/XASW+N6VQWD1hEFwu5S/TYyxewBc0t2LE9HXiIgR0VCxTET0CyJaRURL1NnnykXCIXzxvROMo6dVomLcJcmEo+TSCTdbqBXAlIhM0p0eTFjivahK644c7kbvJWezUsTAojDc3P8isiXXnqPUHNIxvjw1F1Bk4r0sPgd97o9on0N2zUESGKNA/uPjDC6U55BluvV9R+PrFxzp95FoDabeGH5s6iEAgIkj+rvrch3UVZEg7GzuwBW/edU4I56KPhGWaZ6IBd8+D7O0SXnUY+QtnHLoEMz4eoPx3PlAlL1OHlQf1BycCG2mGMTSHLRG2gHXJHKPD/OfcwyA8wGsV1ZfBGCC+DsZwL3if48njl27QmkAEo5nYgpmzvR+R31AcSuKqYnTD2VZzEqApzmY5pjIRjbz0MAsUTa5NPDeoLr4g+BU5JiHdMikKmdMGIpX3tkJIAeHdMhGN1op4Z/7w+9z0KKVIuqaHPAl9zb7HJTlbpiVdL+Ib5vmEL9k8khcMtnfn8y1kVX3f+mtbZGRebrm8NSXTgtogmECmx/j1Tld8OUagmsiqq7UVyWN3z1R9zpmuRL37fxM+fsRgJMAfKSb174bwI3w19/LATzEOHMADCSioI2iBxLnlSUVgZBwwl+02uvKJ6+8zh0fmIxzJw3HMcp0iHrlTGfMA4tUZKPZnXTyYY2kHoIpkbvnOkIa8DSHOFqdinwvYZrDw5/x+iuq5iCvEpU/Kmx9hRM+2Y/uJ4gy/emJE4NCVtMccnBI642kLsDCtoWRayOr7t/Smca+9vDJKP2zJPK6a0r3EoZDXtZk/V5yDcE1EWkSC5lBMiqCqhjEjVYy6155QkSXA9jEGFusVfRRADYoyxvFuh6fKly9jzs/NNm4jzpCWjcVqKgfb2Ui9xwvOkeOqMfvr/LPzRQYIR2RPkMizS3dmdYj7Pzh4wD4/1y+CX1Kyahj77vyRLy9rdlYllg+h2S8aKXwdy3KnNBDWf0+rrihrIEoIT1vF8JDZsPQkzpKfMIhYFaS1w8nV7OoXt+jpsvMNX9U4HiHQs2ohei9Z5v7Q15C36uUmkNcs9JXo7Yzxu4yHDMdfA4InVsAfAvcpJQ3RHQdgOsAYPjw4WhsbMzrPM3NzXkfq/LfV15xf/dvWoXGxtWBfWbOnImOdj43xK4dOzBvbhMAnrpCLcOy7V5WxrdXLAu95huLFmDP6ugEYiaam5sxb+5c37pFi95AS0sHdu5oD30eazfyntqu3btzfmZdXfzYRQsXYOc75kbh+uOqMLKW8MyaLlxyaAqNjY1ob2sHAMyZ8zrW1sZrTNas49eSs9bNmjUTrS0txjJXA5icABobvclpVm7mz3/L1q1Z73PunNluSOY76/l1Tddq6fIaMnXb2jW8jHt27cSbS5rc9UveWOTbv6nDO37L5s1obNzlO/9XTqzC+v0ZrFzxFgBg9+7daB6WQjrtb0yWLFmM9Xu8dCjz583DxizPtbmFp3NYvGgB9ir1bd1ar+e+cMF8bKvzzrNyC3+GHR0doc/w7a3B+SeinvfuXe2+5a6ursD+8nve1uKFe86bOxfrY9YdyYb169HZxp/TmtWr0Zj2rN+vzp6Vs8DZus0r+/79zZg/b17ovqmuLixevBgAsG+fl0yxsbHRbT8kW7ZuRXO/4HMoBLlEK00F8JRYvgzAXADvhB3AGDvXtJ6IjgUwHoDUGkYDWEhE0wBsAqBOEDBarDOd/34A9wPAlClTWENDQ8xb8dPY2Ii8jn3uad/imWeeCUx/DgBw1pln+tV7sW9DQwPqFjRiW2sLRo4YjtNPmwQ0vgQi8pWBrdgOLOSV58TjjwMW+htyyaknT8MRBkdbNhobG3HssVOBVxrddZOPOw7Va5ZixPCBaGg4wXjcvsWbgTcXYcCAQWhoOCWna1bOegEtXV2YNnUqjhxhLnOD+H+Vsq5m3gygrRXTTj4Z4w0pNkysf20doAjVsxsaMHPmzNjvef/izcCSRRgy9CA0NITERIh3eu7ZDa7GsPn19cDypaiprUVDw1n+c7Z3AS+9AIDXg2tbliPDRLqLFctx8IjhOOGEMcD81wEA06aeBLz+X3f/Xc0dwIzpAIDRo0ahoeEY3/nlnf1z0UZg6WIMGTIEdXUtIKfNN0fnCccfj/S63cCqtwEAp55yMsYOiX6u/RY0Ai0tOHnaVJ9DeU1yLbByOQDgPaf430/r0i3A4oWorKwKfe4dy7YCbyzw30fEO3py6yJg62ZvhZMI7C+/5w27W4FXZvCynXqKcbpOlZ/Wb8TX/7bYXR4/biy2s91Y07QbE4+cgIZTx7nv/L1n525I+de2N4DNvCmrr6/DySefAMyeady3X3UlTjrxBGDea+jfvz+wdy8A/mzqF83ElhZP0x0xYgTq6vbk14ZlIa5wGA3gRMbYfgAgolsBPM0YuzLXCzLGlkJJ2kdE6wBMYYztJKKnAHyRiB4Fd0Q3McZ6vEkJiA4TVVFz6bi5lbR9TLM/RZ0rH4zpMzJZHNI5RPGEkau51nWs5nDN4NzCuT2nXMxKvgnuIy6jm0RuueQoAMCDr67j2xNOYFIaU5n0a+oE3p/J55DrGIA8HNLuOIcIw1JU3Tah2/pTEWYln88mxnfyoZNG497GVVi9g4erOkQ5jXfJlagosYQSKaW/a92013DkMGD3noKXD4gvHIYDUIfpdop1heYZABcDWAWgFcCni3CNohPHJhyVlTVutJJpGst8ybDohGBA93wOspLn2lDXu3H78Y/r7sfsCocc7zPMTgyERwXJYyoS5Ev+pwuTuA26+5zFsh5ZpudWihOzH5ajKsrnEM8h7YgyxYvX1x3YcVO8x50LRC2CTzjk4bPQ0dPaZHuHYVulEB43pAZP/+8ZqK1KorHx7W6Xz3itmPs9BGAuEf1TLL8fwIOFKABjbJzym4GnBz/giOsAUx3SceZzMM2MJclnMI57De3accY5uCGe3QhXytVW+9tPnoR/LtqEcUOizQIq3RYOlF1z+MHlR+PV1X67f1RUVFgDIy+RcMjXEQg0tr76FXoZJWLKvJ1HK2Uvl4p8ClHCIZCVVTqkI6qKOld3VEPv7p9DZ8j3vGIeppbVIa++y3tTZ3LMla9fcCSaO1J4dB6Pt8lX3siyEFHoAL5CETda6XYiehbAGWLVpxlji6KO6Sv8+4unY1Bthe9lR8a7G6KVTMnRJPlMVxiHfMY5yLIUI5Q1jIMH9sMNZx9e1GvonD5hKC4+dgRuunBS6D6fPHUcPnnqOP/KGO9dR6agrkg4muZAaPx6AzbuaQscH6Xdee9PmnW0ImqaQxyTS9ikSb65BiI0nTB80+6GzxDrkkt91wfBxUHVshxFWEtt7K/XnsJ9R3lQXZHA/5x8iCscTM8n4ZA7w13YJ5Z0hUNexciJXERPDYB9jLE/EdEwIhrPGFtbrIIdKBw7mo8dSCk9nzg24ajcLf5Q1uKYlfQPRo5ziC47/5+XWUn+L0Gl7q7mUF2RwG8+cVLOx8lnarrHsOcqw20rEo6vI5BwCKMH1WCccPLG7XyYhL6+XT0+3ghpcV1N+qnH6nN8ez6HcKSZKG54Zi6hr2p54nYWfMJBMStJAd6vMuEb15Ir6vMz3bIqHMJwNYe8SxGfuKGs/wcesXQkgD8BqADwCPjscBbEt4nLD0HONDV2SA2+fM4E3z4+4RA1QrobmoN+ZIYxZLKkJJDl6p5DuvjVupQDhVTyuWpKmYNZdUjrDaHfrBRlr/Yv61qePmYiTh1yhYO2q9qT1zsq3gjp8PPKe5TnyfYJ5RLjn4ypaakoQV1wyPv2OrpiqDU5wJhZYCUdQieizZPyWeUzdiNX4oriKwC8D0ALADDGNgPIPYayFxO33qrz4RIRZn7j7EBqYvVcUb2l7mgOujCTZqU4ycwyeWQMlpcrRaUum3DI47JyIFfSoDmoxB0EJ0VUtM8ht171dy49CkNqK3GQlkZaNd2E+Ryi8GZWdJSSh5PLyORCTPYjNYc4/pBcMY6iJ9nwK/tp+0Rpp4UmrlmpkzHGiIgBABHFCzjvQ8TVHNSsrGH4QllL5HPg04TGS+o2cWT+/YJStNvlEg7uh5uDDiHnJ6hw9Ggl/zlIhDeGzfOdSxl9DukYz+rCY0bgwmOC41nVWfDC03mHqw6y41Phag7RZcnlrhM5CkDAb1YiIkVzKKxwIDI37okYz8H1OZTAsBRXODxORL8Fz3V0LYBr0Acn/smXQTUV2CNmc3Md0jFj1SOjlQo4zuHx+Ruyps8YM7gGj113CiaPHpj3dUthLC1E6GE+5HPZlKI5+CfPCb73BBFSWeYxcMsSUUb1+O4IGnVuZr2zE0foJDWzUrZDchrr4jMrxTtGFQ4J8qKVOlKlmVzHa/jDcbWsnqA5EK9JjwGYCGAfuN/hu4yxF4tctl7Dazef4/5WzUph+MY5RKjSuY4Z8B/rX25cuUNcO/qcJx86JO9rlopyaQ75IM1KFQkt8Z7hHvi7ifYLZSOf3EphRDln44Syeg5pf1RQGPl6uuJHKynHOOTOiDhmsHkSre5gei5uOSn8uXnRSj1AcxDmpGcYY8cCsAIhD9RUGvKjiReOGN+Zlithlat47apsLYp1fo+ym5VyuLw0K0WNPnaJ5beJTlPOp6dVl/N/VlFzjcRpkKXPLG4ETj6p4oEcBsFpZqWLjx2Bxz93KqaOG5TXdeNeSxLnOSRialmFIK6HZyERTS1qSfoIYbnxVUqhMoZdvlgO41Jaeg5Uh7SKqX7INVG3FxZ2KuGT2EhHcPeeU7+K7MIhMpQ16c96mu355RskF1cA6oPgiAjTxg8uSi99cG0lJmo5xvwD3PizPVib+reU4xziCoeTAbxGRKvFDG1LiWhJMQvW24mrORTt+mGaQ5Ea1m+cfyQAbzKaYlKK52ciHyfhSWN5r3SS5uQ3NUieZhJ+nWztJynn6bZwiNAcwgZ4qsiGLm54ZncGX8YhX80kH5IJB8995Uz/OkVzOPrgAfjlx0/ADz9wrG8fqQWV3SGtDHS7oOgl6SPICtAZ4eSK89GeO6l7qa2iQh2LwUemjsFHpo7JvmMB6E4UV3fI59l98MRROO3wIRg5ILtd2zVJxnFIh+xCSrRSdx333TYruSOk44WydmcekTiUTjSYkR0z+ewuO+7giH1KUJ4s2/8u/v+RMfau/lfswvVGPnvGoQAQmX46zovXJ+/JlVDNoUy97kJSNs0hj8sSUSzBIPcFutcwOORP4dId1ImOAtcRLUtUg5twuKCq0BrFMIrdeGcU1aQYckjenp5ZVRLHZOQK9B7gkHaI6FsAjjBN+GOa5McSzamHDcHK2y70jYbV0c0GP/ngZNz4j8Ja8aJ6lgc6pZwty09xrxtnIGG2Rk31OXT3MUXOwRyzHlUkvKlzsz2+TJHtSuqzK8aVjhrZH9eeMR6f0nNyCcKEhorcpRRVPJtw+Bh4BtYk7IhoH7/4+AnY09KZfUcDUYIBCH78H5k6Bn+YvRYrt+3P63pxruGtL9glykb5opWKff74zshw4e+VsxBzIYdfX/ocoverTDhe6u4s59RlQ9zJn+JSbJ+D45A7j4cJ+Tr0DtqzXz4D+9r4OClvoGXxySYcBgP4CYAljLFnS1CeA4b3GeyBhcIYxVjg2hB2unINICsk5TMrZXcYdwe3gx3pkI5u4PgI6XhmnO4QxyENiLmzlXxjUaiNd//qJJ798hkRe+eOT/iU0DktSYQ0/JNGerPvqRFNxSZb1+FTABYAuIqIriYi05zQlgJj+mgL/SGXOlqplJTLIV3sqzo5mINkNMsfrpqCEw4Z6K0nNVqp4EV0ycWsJMlWx7/03sMxVszrUVeV9E/FWwBUwVoO53Qihs8hlzrQXSKrB2Ps84yxEwHcCmAQgAeI6DUi+iERnUlEhX07FgDmylHASd9CrxG1/kCiJzukZX6q/M6fvcevd3jPmTQcv73SSz+u5lYqppYoT52tka2pTLhZALKZA4fUVeFekUq9GD3nYofKZiMRwzGfT/6ufIk72c8KACsA3E1E/QCcDeDDAO4CT+VtKSCmj6TQDZ76cakzXPWGaKVyOaSz2YOXfe+CrM/3EycfgoXr94ac3//fhNu+Kfvo81y7mkMRNay4fp+ff+wEDK6txGEH1eKDWnZiE/pUoQVFdUiXQVDE0Rw8v0TxyxMnt1ICwDLG2EQAYIy1gc/1/EyRy9ZnMTUgxbQxPnnDafjQva9i/rt7rEO6O2S5bJxpHW+/4tjQbV4UY/b7U/fwp/smVwsthObwm0+c6E6Go+LEVB2OGzMQAHDzReGz7qkU04muzlNS7DEVkulfPQvn3jUTgJJjKmJ/J4YAKRRxciuliWglER3CGFtf/CJZTC++6PbsmLHmBwLl8puUzucQZVYKNmq++QFIMU8V4DldfOxI4/o44xzyoZhaoerwLpXicPhBde5vL6T3ADIrgfsblhHRXIgJfwCAMfa+opSqj2N2SBf3msleJBzKblYq0uVzckj7TEl+zUFqDMX0ORyIUW/HHDwASzc1ASizWSlqH1ewF788cYXDd4paCouPUkQr6XjOsKJepiT0ZId0Ic6f6/2pZiVSji+m+c1NvFeOVjZPHv7MNFz5h9fx5qZ9Zbl+rGilEk72E0v+MMZmAlgHoEL8ngdgYRHL1acxfbPFbvCcApoayk355pAuzTvKtSqoj8MXrVRM4VCkc/ev5tFe5046qODnHlhTiSljBwMoUyhrjAFupazasTQHMfvbdeCD4g4DMArAfQDOiTrOkh+mD6vYvdJSTiJSbMo9Qrp4ZiX5PzeHtKOZlQqVsjsKd7KfAp93QE0F5tx8DobWVRb4zBzXj16OQXAxTLtSgGQb7FgI4lqubgBwGvhMcGCMvQOg8KLbAqA8ZqVSZnssNj01Wqm7ODHszZEzjAEAxWuEuotnVir8uUcMqC5a1FIpzDVh5GJWypRg5tK4T7iDMeYmEiKiJMqf4bbXYjQrFdkBVQonZakol0O66A1LDM3hwmNG4IoTRuGWS7zQUL9ZScmtVMTn1BOqUf/quC7VnkEyhj/BKaHmEPfpzRTZWfsR0XkAvgDg38UrVt/GOM6hyA1PKXqTpaJc9+AlFy3O9eNM9lNdkcDdHz0eAPCWdpz8XchQ1jBkzP65R3Vv3pHuMOdb5yCV47Bnz6xUhAJlwfEqUChSYSpF+eL2R28CsAPAUgCfAx8A9+1iFaqvYxznUGyTRQkH1xSbsmkOxY4oyzOvjioESpU+ozLp4NWb3ouffnhy0a6RjZrKpOvAjsvUcdwhfcyoAcUoUiTqTHBhFNNcFyhPzP3OBvAIY+x3xSyMhVMOn0Nvckg7DuGPV0/BNQ/ML+l1e2ooq36OQk32kw19/uMDgQuPGYH53z4XQ+uqSn7teOkzep5D+lMAFhPRHCK6k4guI6JBxSxYX6Ycg+AOxLj0KN47sfTmjFK9o+5cx5+V9cDvCBSDcggGQB0EF8Pn0FM0B8bYVQBARAcD+BCAXwM4OO7xltwIG+eQcAjpDMOI/tUFv6a0ZZZykvXeR3FNc15QQv4X4D4H/tsKh55FnGg0N1qpBN9p3HEOVwI4A8CxAHYC+BWAV4pYrj6NybRDQjgs+94FBWt87v7ocTh0KM/tIhuKXB14Fo+i+4UKoDmUarIfS+7E0Ry8cQ7FJ27P/x4Aq8EHvs1gjK0rVoEsZhzifoFCTnByxQleimTZUBR7nt7eTLGbWnn+bvkclOPLN9e2xUScQZTu4MKeEq3EGBsK4BoA1QBuJ6K5RPRwUUtm8SHNSsUiaTWHbuNOE1rk83dnzAt3SPPf1qzUs4gTDOI4pdMcYlUzIuoP4BAAYwGMAzAAQAnG6FkkapRJMZCVLm2FQ48lzhzS2aASjXOw5I58G1Hv1zUr9RSfA4DZyt+vGGMbi1ckiwmHqKhmAFnpepND+uPTxqA+xzj37lDsD7ZQvgI3WsnKhp5F9jFwXihrD4pWmgwARFSXbd84ENGtAK4FH1gHAN9ijD0jtt0M4DMA0gD+lzH2fCGueaBTbM0h4WoORbtEyfnRB8o3AKsYFGpyeS8rawkmBbDERjqio96vN8FeD9EciOgYAA+DZ2UlItoB4CrG2JvduPbdjLGfatc5CsDHABwNHio7nYiOYIylu3GdXoE6SUtRzl/CEDlLfhRiEJx6fBFn3LR0g0izUg9MvHc/gK8yxsYyxg4B8DWxrtBcDuBRxlgHY2wtgFUAphXhOgccSYdQkSze1+w6pNNWOOSL++SKJMQLFXoqOwLWId0ziTQr9cBQ1lrG2Ay5wBhrJKLabl77i0T0KQDzAXyNMbYHfJ6IOco+G8W6AER0HfgcExg+fDgaGxvzKkRzc3Pex4ZRyPPJcx3fL41x41hBzm265/Xv8qS7a9auRWPjpm5fo6dRjPess2ovV3D379tXlGs1NbUBAN54YzG6NmYPaTbdc2NjI3a3827n1i1b0Ni4u+DlLCeleM/FYsP69QCAXbt2hd7D8m0pAP77LNY9xxUOa4joO+CmJQC4EsCaqAOIaDqAEYZNtwC4F8APwAXgDwD8DDxUNjaMsfshtJcpU6awhoaGXA53aWxsRL7HBnjuaQAozPkKeS4N0z0vY6uA1SsxduxYNDQcWfBrlpuCvucQ+q/fA8x5FfX9+6Oh4bSCn//3q14Hdu3EsZMn48wjhmXd33fPSn3atq8daHwJh4wehYaGYwpeznJSivdccMS7GTv2EGDtagwdOhQNDVOMu3Yt3wYsmo+amlo0NJwJoHj3HFc4XAPgewCeAG/QX0GWxpwxdm6cExPR7wD8RyxuAjBG2TxarOuTjBlcusRlV79nHDbuacO1Zx5asmv2NqS7pnjjHMR1unme3jQlbG9Cvt/oxHv8f9kd0kRUDeB6AIeDp+v+GmOsq7sXJaKRjLEtYvEKANKx/RSAvxDRXeAO6QkA5nb3eqXigU9PRf9+hQmdfOXGswt2rjjUViXxow8cW7LrWXLHHcXezaCBUqTstuSOjFaKeivSWd0TQlkfBNAFrilcBGASgK8U4Lo/IaLjwTtB68DniABjbBkRPQ5gOYAUgBsOpEilhiMLN3PqmME1BTuXpXfgpU7ornCwDumeSJxotJ7kkD6KMXYsABDRH1CgXjxj7JMR224HcHshrmOxlIOiZWV1818V5jxWOPRM4sznUIqQ82yxka4JiTGWKnJZLBZLBFQgsxLZ3Eo9Ei99Rvg+TqEcTzHIpjkcR0T7xG8Cn0N6n/jNGGP9i1o6i+WAotjpM/j/7goHb7pRKxx6ItGT/fD/ZTcrMcYKlx/aYukjFKvJ9RzShTmP1Rx6FucfPQKrdjTjW5dMCt2nlGYlO5ObxXKAIFMhdbddSDiEyqSDuir7+fck+lUm8JtPnBS5j1PCxHs2u4rFcoBw7KiBAIARA7o3TWxl0sETn38PPjp1TPadLSUjjpnPnUO63OMcLBZLfIrdm/vcmYfijAlDccyoAd0+VyHOYSkscax8iQJFrMXBag4WS4E4/CCe0f7aM4ozytxxyDbqvZhYmkMJW2yrOVgsBWJgTSXW3XFJuYthOUCJEzzmTfZT/nEOFovFYikBseaQLuEIaSscLBaLpQcQx+dQqPxacbDCwWKxWHoAcXwOpZxD2goHi8Vi6QHE8TmUcoS0FQ4Wi8XSA8hpnIM1K1ksFkvfwJqVLBaLxRIgF4e0NStZLBZLHyFWKKtosW20ksViKQiTRw9AdYX93HsysQbB9aBpQi0WSy/gqS+eXu4iWLKQm8/Bag4Wi8XSJ4jjcyDrc7BYLJa+RbxQVv7fRitZLBZLHyGOz4HsOAeLxWLpW8TRHOQeZZ9D2mKxWCylIY5w6FeRwISD6vCVc48oenmscLBYLJYeQKxBcA7hxa+eVfzCwJqVLBaLpUcQZxBcKbHCwWKxWCwBrHCwWCwWSwArHCwWi8USwAoHi8VisQSwwsFisVgsAaxwsFgsFksAKxwsFovFEsAKB4vFYrEEsMLBYrFYLAGscLBYLBZLgLIJByL6EhGtIKJlRPQTZf3NRLSKiFYS0QXlKp/FYrH0ZcqSeI+IzgZwOYDjGGMdRHSQWH8UgI8BOBrAwQCmE9ERjLF0OcppsVgsfZVyaQ6fB3AHY6wDABhj28X6ywE8yhjrYIytBbAKwLQyldFisVj6LOVK2X0EgDOI6HYA7QC+zhibB2AUgDnKfhvFugBEdB2A6wBg+PDhaGxszKsgzc3NeR97oGLvuW9g7/nAoqe1YUUTDkQ0HcAIw6ZbxHUHAzgFwFQAjxPRobmcnzF2P4D7AWDKlCmsoaEhr3I2NjYi32MPVOw99w3sPR8Y/GP8bqza3oyGqYfkdXyx7rlowoExdm7YNiL6PIAnGJ8IdS4RZQAMBbAJwBhl19FincVisfRKTho7GCeNHVzuYgQol8/hSQBnAwARHQGgEsBOAE8B+BgRVRHReAATAMwtUxktFoulz1Iun8MfAfyRiN4E0AngKqFFLCOixwEsB5ACcIONVLJYLJbSUxbhwBjrBHBlyLbbAdxe2hJZLBaLRcWOkLZYLBZLACscLBaLxRLACgeLxWKxBLDCwWKxWCwBrHCwWCwWSwDiEaQHNkS0A8C7eR4+FHyMRV/C3nPfwN5z36A79zyWMTbMtKFXCIfuQETzGWNTyl2OUmLvuW9g77lvUKx7tmYli8VisQSwwsFisVgsAaxwEJld+xj2nvsG9p77BkW55z7vc7BYLBZLEKs5WCwWiyWAFQ4Wi8ViCdBnhQMRXUhEK4loFRHdVO7yFAoi+iMRbRfp0OW6wUT0IhG9I/4PEuuJiH4hnsESIjqxfCXPHyIaQ0QziGg5ES0joi+L9b32vomomojmEtFicc/fE+vHE9Hr4t4eI6JKsb5KLK8S28eV9Qa6AREliGgREf1HLPfqeyaidUS0lIjeIKL5Yl3R63afFA5ElADwawAXATgKwMeJ6KjylqpgPADgQm3dTQBeYoxNAPCSWAb4/U8Qf9cBuLdEZSw0KQBfY4wdBT717A3iffbm++4A8F7G2HEAjgdwIRGdAuDHAO5mjB0OYA+Az4j9PwNgj1h/t9jvQOXLAN5SlvvCPZ/NGDteGc9Q/LrNGOtzfwBOBfC8snwzgJvLXa4C3t84AG8qyysBjBS/RwJYKX7/FsDHTfsdyH8A/gXgvL5y3wBqACwEcDL4SNmkWO/WcwDPAzhV/E6K/ajcZc/jXkeLxvC9AP4DgPrAPa8DMFRbV/S63Sc1BwCjAGxQljeKdb2V4YyxLeL3VgDDxe9e9xyE6eAEAK+jl9+3MK+8AWA7gBcBrAawlzGWEruo9+Xes9jeBGBISQtcGO4BcCOAjFgegt5/zwzAC0S0gIiuE+uKXrfLNU2opUwwxhgR9cr4ZSKqA/APAF9hjO0jIndbb7xvxqfQPZ6IBgL4J4CJ5S1RcSGiSwFsZ4wtIKKGMhenlJzOGNtERAcBeJGIVqgbi1W3+6rmsAnAGGV5tFjXW9lGRCMBQPzfLtb3mudARBXgguHPjLEnxOpef98AwBjbC2AGuEllIBHJTp96X+49i+0DAOwqbUm7zWkA3kdE6wA8Cm5a+jl69z2DMbZJ/N8O3gmYhhLU7b4qHOYBmCCiHCoBfAzAU2UuUzF5CsBV4vdV4DZ5uf5TIsLhFABNiqp6wEBcRfgDgLcYY3cpm3rtfRPRMKExgIj6gftY3gIXEh8Su+n3LJ/FhwC8zIRR+kCBMXYzY2w0Y2wc+Df7MmPsE+jF90xEtURUL38DOB/AmyhF3S63s6WMTp6LAbwNbqe9pdzlKeB9/RXAFgBd4PbGz4DbWV8C8A6A6QAGi30JPGprNYClAKaUu/x53vPp4HbZJQDeEH8X9+b7BjAZwCJxz28C+K5YfyiAuQBWAfgbgCqxvlosrxLbDy33PXTz/hsA/Ke337O4t8Xib5lsq0pRt236DIvFYrEE6KtmJYvFYrFEYIWDxWKxWAJY4WCxWCyWAFY4WCwWiyWAFQ4Wi8ViCWCFg6VXQkRpkcVS/o0rd5m6AxG9n4i+K37fSkRf17avI6KhEcdPl5k7LZY42PQZlt5KG2PseNMGMWiOGGMZ0/Yeyo0A3teN4x8G8AUAtxemOJbejtUcLH0CIhpHfP6Oh8AHjY0hom8Q0TyR9/57yr63ENHbRDSbiP4qe+lE1EhEU8TvoSKNg0yAd6dyrs+J9Q3imL8T0Qoi+rMQTCCiqUT0KvH5GOYSUT0RzSKi45VyzCai44joCAAdjLGdMe7zekVbWktEM8SmpwB8vACP0tJHsJqDpbfST2QsBYC1AP4feI77qxhjc4jofLE8DXxU6VNEdCaAFvDUDMeDfx8LASzIcq3PgKcpmEpEVQD+S0QviG0nADgawGYA/wVwGhHNBfAYgI8yxuYRUX8AbeApQK4G8BUhEKoZY4uJ6NOiHCr/j4iuVJYPBgDG2H0A7hO5pl4GcJdYv4f45DdDGGMHXH4hS+mxwsHSW/GZlYTP4V3G2Byx6nzxt0gs14ELi3oA/2SMtYrj4uTcOh/AZCKS+X0GiHN1ApjLGNsozvUG+FwbTQC2MMbmAQBjbJ/Y/jcA3yGibwC4BnziJoDn69+hXfNuxthPlftbp23/OXguoX8r67aDCxErHCxZscLB0pdoUX4TgB8xxn6r7kBEX4k4PgXPFFutnetLjLHntXM1gM/YJkkj4ptjjLUS0YsALgfwEQAniU1t4AInFkR0NYCxAL6obaoW57JYsmJ9Dpa+yvMAriE+BwSIaJTIlz8LwPuJqJ/IhnmZcsw6eA32h7RzfV6YckBER4gMmmGsBDCSiKaK/euVlNO/B/ALAPMYY3vEurcAHB7npojoJABfB3Cl6nAXvo4R4h4slqxYzcHSJ2GMvUBEkwC8JnzEzeAN6kIiegw8C+Z28PTukp8CeJz4bFxPK+t/D24uWiga4R0A3h9x7U4i+iiAX4p0220AzgXQzPhENvsA/Ek5ZBaAnxERseyZMr8IYDCAGeK+5jPGPgsu1OYwb8Y0iyUSm5XVYomAiG4Fb7R/mm3fAl3vYACNACZqPf+fA/g3Y2x6nuf9OYCnGGMvFaSgll6PNStZLD0EIvoU+NzXtxjGYPwQQE03Tv+mFQyWXLCag8VisVgCWM3BYrFYLAGscLBYLBZLACscLBaLxRLACgeLxWKxBLDCwWKxWCwB/j9hzeJgvrSokgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Pxx=10*np.log10(pxx[1:,]) # 将单位转换为dB\n",
    "f=F[1:,]\n",
    "plt.plot(f, Pxx)\n",
    "plt.title('Periodogram Power Spectral Dendity Estimate')\n",
    "plt.xlabel('Frequency(Hz)')\n",
    "plt.ylabel('Power/frequency(dB/Hz)')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a64d3cc",
   "metadata": {},
   "source": [
    "找到周期图的最大值。Fisher 的g统计量是最大周期图值与所有周期图值之和的比率。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5001cb37",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fisher_g = 0.03797947482896527\n",
      "100.0\n"
     ]
    }
   ],
   "source": [
    "# 找出最大值所对应的频率\n",
    "pxx=pxx.tolist()\n",
    "F=F.tolist()\n",
    "fisher_g=(max(pxx)/sum(pxx))/2\n",
    "print('fisher_g =',fisher_g )\n",
    "print(F[pxx.index(max(pxx))])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e68529e1",
   "metadata": {},
   "source": [
    "确定pvalFisher 的g统计量的显着性水平,在计算二项式系数时，使用 gamma 函数的对数来避免溢出。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "618c7e62",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.032144002118851e-06\n"
     ]
    }
   ],
   "source": [
    "import scipy.special as sc\n",
    "N=len(Pxx)\n",
    "nn=np.arange(1,np.floor(1/fisher_g)) # np.floor向下取整函数\n",
    "I = ((-1)**(nn-1)*np.exp(sc.gammaln(N+1)-sc.gammaln(nn+1)-\n",
    "       sc.gammaln(N-nn+1))*(1-nn*fisher_g)**(N-1))\n",
    "I=np.array(I,dtype=np.float32)\n",
    "pval = sum(I)\n",
    "print(pval)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
